Charge ordering in quarter-filled ladder systems coupled to the lattice 
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We investigate charge ordering in the presence of electron-phonon coupling for quarter-filled ladder 
systems by using exact diagonalization. As an example we consider NaV205 using model parameters 
obtained from first-principles band-structure calculations. The relevant Holstein coupling to the 
lattice considerably reduces the critical value of the nearest-neighbor Coulomb repulsion at which 
formation of the zig-zag charge-ordered state occurs, which is then accompanied by a static lattice 
distortion. Energy and length of a kink-like excitation on the background of the distorted lattice 
are calculated. Spin and charge spectra on ladders with and without static distortion are obtained, 
and the charge gap and the effective spin-spin exchange parameter J are extracted. J agrees 
well with experimental results. Analysis of the dynamical Holstein model, restricted to a small 
number of phonons, shows that low frequency lattice vibrations have a strong influence on the 
charge ordering, particularly in the vicinity of the phase transition point. By investigating the 
charge order parameter we conclude that phonons produce dynamical zig-zag lattice distortions. 
A model with only static distortions gives a good description of the system well away from the 
transition point while overestimating the amount of charge ordering in the vicinity of the phase 
transition. 

PACS numbers: 71.10.Fd,71.38.-k,63.22.-|-m 



I. INTRODUCTION 



The formation of an ordered pattern of ion charges 
is a rather general type of phase transition which oc- 
curs in three- as well as lower-dimensional solids. It has 
been known for more than six decadesi since its discov- 
ery in magnetite Fe304. Even in that compound this 
phenomenon still attracts a lot of attention due to the 
interesting physics of the transition. Since the charge 
ordering causes changes in the interaction between the 
ions, it drives a lattice distortion, which, in turn, in- 
fluences the ordering pattern. The quarter-filled ladder 
compound NaVzOg (Ref.Q) is another interesting exam- 
ple of a system which shows charge ordering. Here the 
transition, as observed in the nuclear magnetic resonance 
experiments;^ occurs at Tqo ~ 35 K and is accompanied 
by the formation of a spin gap^ at the samc'"^ or slightly 
lower temperature. In NaV2 05, where one d^y electron 
is shared by two V ions in a V-O-V rung, the ordering 
occurs as a static charge disproportion S between the V 
ions, which obtain charges 4.5 ± 5 with a zig-zag pat- 
tern of (5's. Most probably, the main driving force for the 
transition is the Coulomb repulsion of electrons on the 
nearest-neighbor sites within one leg of the ladder. For 
half the ladders the apical oxygen ion is located below 
the ladder and for the other half above. Therefore the 
crystal environment of the V ions is asymmetric, and the 
dxy electron is subject to a strong Holstein-like electron- 
phonon coupling.^ As a result, the transition is accom- 
panied by the displacement of ions from their positions 
in the high-temperature phase (T > Tqo)- They then 
form a larger unit cell, as clearly seen directly in the x- 
ray diffraction experiments^ where displacements of V 
ions of the order of 0.05 A were observed, and indirectly 
in the appearance of new phonon modes in the infrared 



absorption^ and Raman scattering spcctraiS 

The importance of the coupling to the lattice for 
phase tr ansitio ns in quarter-filled systems was shown 
in Refs. llOilll The low-energy excitations of the zig- 
zag order parameter are also strongly influenced by the 
latticeJ^ 

The static properties of the ground state of the doped 
ladders without coupling to lattice distortions were ex- 
tensively investigated using mean field approaches 
the density- matrix renormalization-group (DMRG)^ 
bosonization, and renormalization group techniques^ 
The DMRG studies show that at strong enough Hub- 
bard interaction the ladders exhibit a zig-zag charge or- 
der if the repulsion between the electrons on the nearest- 
neighbor sites exceeds some critical value Vc- The corre- 
sponding phase transition is of second-order as a function 
of V. For a very strong intersite repulsion, phase separa- 
tion becomes possible. Dynamical properties were stud- 
ied with exact numerical diagonalizationi^iiSiSS without 
taking into account the coupling to the lattice. They 
were quite successful in understanding NaV2 05 dynami- 
cal properties above the transition temperature. 

The investigation of the ordering of electrons interact- 
ing with the lattice requires the knowledge of electron- 
phonon couplings and lattice force constants in addition 
to electronic parameters such as hopping matrix elements 
and electron correlations. For a given compound al- 
most all of these parameters can be extracted from first- 
principles band-structure calculations. The force con- 
stants and electron-phonon coupling can be obtained by 
comparing the total energies and the intcrionic forces in 
distorted and undistorted lattices. The phonon frequen- 
cies required for studies of dynamical lattice distortions 
are given by experiment^ while the necessary knowledge 
of their eigenvectors can be obtained from first-principles 
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calculations. We performed such calculations of the band 
structure, lattice dynamics and electron-phonon coupling 
for the NaV205 compound. These calculations are in 
good agreement with Ref. 2 in the part concerning the 
band structure, and their details as well as a comparison 
to other first-principles calculations^^ will be published 
elsewherej^S In the present paper we will concentrate on 
the strongest electron-phonon mode present in NaV205, 
which is a simple Holstein-type interactioniSS 

We therefore investigate a model which takes into ac- 
count the main interactions, i.e., the Hubbard and inter- 
site repulsions and the coupling to the lattice. We study 
the ground-state properties of a quarter-filled ladder cou- 
pled to static lattice distortion with the Lanczos algo- 
rithm. We show that the ordering is strongly enhanced 
by the interaction with the lattice for realistic values of 
the coupling. Secondly, the interaction with dynamical 
(quantum) phonons at different frequencies within the 
Holstein model is considered. We find that indeed, the 
quantum phonons produce and stabilize the static order. 

The paper is organized as follows. In Sec. ^1 we in- 
troduce the extended Hubbard model (EHM) with addi- 
tional lattice distortions. In Sec. IIIII we present results 
for static properties of this model, including kink exci- 
tations. Dynamical quantities are discussed in Sec. IIVI 
Section^focuses on the influence of dynamical phonons, 
and finally we give our conclusions in Sec. IVII 

II. MODEL 

The quarter-filled ladder compound Q;'-NaV205 can be 
described microscopically by an extended Hubbard (or t- 
U-V) model (EHM). For the description of the distorted 
low-temperature phase we also include the coupling of 
electrons to the lattice, yielding the model 

H ~ Hehm + Hi + He-i, (1) 

where Hi is the lattice deformation contribution, and 
He-i the electron-lattice interaction. These terms are 
given by 

+ U ^rii-friii -I- ^ VijUiUj, (2a) 
Hi=nJ2Y' (2b) 

i 

He-l = - C'^ZiUi, (2c) 
i 

with the effective lattice force constant k and the Holstein 
constant C. The sites are labeled by the indices i, j and 
Zi is the distortion on site i. The hopping matrix elements 
tij connect nearest neighbor sites (ij) (see Fig. ^ with 
occupation numbers Ui = n.ii + . 



The first-principles calculations done in Ref. |3 and by 
our group^^ yield the intrarung hopping ta ~ 0.35 eV, 
which we will use below as the unit of energy. For the 
hopping along the ladder we use th = ta/2, again in agree- 
ment with the band-structure results f^iS^ whereas previ- 
ous DMRG studiesi^ were mostly done at ta/tb < 1.4. 
For the on-site Coulomb interaction we use U = 8.0 
as estimated in Ref. .2,, We assume Vij = V to be the 
same for all bonds, and take as a free parameter since 
there is no unique procedure of extracting it from the 
band-structure calculations. The lattice distortions are 
expressed in units of 0.05 A since the ion displacement be- 
low Tco are of this order of magnitudeii With the chosen 
units of energy and length the comparison of the band 
structure and lattice force calculations done on distorted 
and undistorted lattices give the dimensionless constants 
K = 0.125 and C — 0.35, respectively.^^ The effective 
coupling parameter C'^ / k is close to unity and, therefore, 
the lattice plays an important role in determining the 
properties of NaV205. 

The Hamiltonian in Eq. will be used for calcula- 
tions with both static and dynamical lattice distortions. 
All quantities presented in this paper are calculated by 
the ground state Lanczos method on single ladders of up 
to eight rungs, which enabled us to perform simple fi- 
nite size scaling. The largest Hilbert space for the eight 
rung lattice considered in this study was of dimension 
-^states = 1656 592, which could be reduced in special 
cases by exploiting translational invariance and 5*2 con- 
servation to A'statcs = 103 820. The next lattice size ad- 
mitting charge order would consist of ten rungs, which is 
far beyond our computational capabilities. The restric- 
tion to a single ladder is a considerable simplification 
compared to the structure of NaV205, but for the quan- 
tities of interest in this paper the role of other ladders is 
of minor importance, since the frustration of the ladder- 
ladder interactions significantly reduces their effect in the 
zig-zag ordered state. 

III. STATIC PROPERTIES 

A. Charge order 

To investigate the connection between the lattice dis- 
tortion and charge ordering we calculated the static 
charge structure factor 

= ^ E e^''^^'-^^-) iin^n,) - (n)') , (3) 

where N is the total number of sites in the system. The 
zigzag charge order parameter mco can be expressed in 
terms of this structure factor as 

"^co = ^^^c(Q), Q = (7r,^). (4) 

The term (n)^ in the denominator ensures that the order 
parameter is equal to unity for full ordering, which in 
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FIG. 1: Schematic picture of ladders, with hopping matrix 
elements ta on the rungs and tb along the chains. The dark- 
ness of the circles corresponds to the charges on the sites of 
the ladder. The upper ladder is shown without charge order, 
the lower ladder with zigzag charge order. 



NaV205 corresponds to the charges of V ions within a 
rung to be equal to +5 and +4, respectively. 

At this point the lattice distortions Zi in the Hamilto- 
nian are external parameters of the model and not dy- 
namical variables. Therefore they have to be fixed in a 
proper way, for which we chose a mean field approach. 
Considering the distortions as the mean field parameters 
one can extract the optimal value for the Zi by look- 
ing for the minimum in the ground-state energy with re- 
spect to Zi. This procedure could be done within the 
unrestricted Hartree-Fock approximation, but this com- 
plicates the calculation because of the larger number of 
variables for which the minimum has to be found. In- 
stead we restrict ourselves to a single order pattern, the 
zig-zag oiAei,^^^ which was observed experimentally^ 



ze 



(5) 



and investigate the total energy as a function of z. 
The optimal values of z where the total energy reaches 
the minimum for several values of the nearest-neighbor 
Coulomb interaction V determined in this way are indi- 
cated by arrows in Fig. El In the following we denote the 
position of the minimum in the ground-state energy by 

For interaction strengths of = 1.5 up to F = 3.0 a 
clear minimum occurs at z « 1. Additionally one finds 
a maximum at z = 0, which results from the z — > — z 
symmetry of the system. We also found a small distortion 
for V = 1.0, but it is strongly size dependent and rapidly 
decreases with increasing length of the ladder, whereas 
the distortions marked by arrows in Fig. |21 are almost 
independent of the system size. Therefore we argue that 
the finite value of z at F = 1.0 is due to finite-size effects 
and disappears in the thermodynamic limit. For this 
reason we did not mark it with an arrow in Fig. |21 




FIG. 2: Ground-state energy per site as a function of the dis- 
tortion 2 calculated on an 8 x 2 system with periodic boundary 
conditions along the ladder, and C — 0.35. From bottom to 
top: V = 1.0, 1.5, 2.0, 2.5, 3.0. The arrows indicate the 
position of the minimum. 




FIG. 3: Charge order parameter for several values of V. Up- 
per/lower panel: Calculation without/with coupling to the 
lattice. The solid lines are obtained by 1/A'^rungs finite size 
extrapolation. 



This behavior gives a first idea about the charge or- 
dering of the system with and without coupling to the 
lattice. From Fig.|21one could expect that the charge or- 
dering transition in the presence of static mean-field-like 
lattice distortions occurs in the region between = 1.0 
and = 1.5. In order to investigate this transition we 
calculate the order parameter given by Eq. This 
quantity shows strong finite-size effects, which makes it 
necessary to apply finite-size scaling. We calculated the 
order parameter for systems of four and eight rungs, re- 
spectively, the largest system size available. Although 
higher-order corrections to the scaling behavior are ex- 
pected for these small systems, we performed a 1/iVrungs 
extrapolation to l/iVjungs = 0. This procedure does not 
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FIG. 4: Optimal distortion Zopt as a function of the Holstein 
constant C for different values of the coulomb interaction V , 
calculated on a 8 x 2 cluster. The horizontal line indicates 
the experimental result. 
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FIG. 5: Contribution of the lattice energy Eq. lllbj (dashed) 
and the electron-lattice energy Eq. I)2c|l (solid) to the ground 
state energy as a function of the distortion z, &t C = 0.35. 
The lattice energy is independent of V . Electron-lattice en- 
ergy, from top to bottom: V = 1.0, 1.5, 2.0, 2.5, 3.0. The 
arrows are drawn where the total energy has its minimum, 
as in Fig. |2l 



give the exact value of the order parameter in the ther- 
modynamic limit and does not allow to extract an exact 
value for the critical Coulomb interaction Vc, but it pro- 
vides the possibility to obtain a rough estimate of the 
interaction V at which the phase transition occurs as 
well as the approximate mQQ(V^)-dependence. 

To study the effects of the lattice distortions we also 
calculated the order parameter without coupling to the 
lattice, that is for C — 0. The results are shown in Fig.O 
As one can see in the upper panel where no lattice dis- 
tortions are present, the order parameter changes rather 
smoothly when going from the disordered phase into the 
ordered one. A different behavior can be found for finite 
lattice distortions. Here the charge ordering sets in, at 
a much lower value of V than in the absence of distor- 



tions. In addition the transition sharpens considerably. 
For both zero and finite distortions the finite size scaled 
order parameter is slightly negative for small interactions 
V , which is due to the fact that higher order corrections 
in the scaling have been neglected. 

The optimal distortions z for different values of the 
Holstein constant C and repulsion V are presented in 
Fig. ^ In our units the distortion found experimentally 
in the charge-ordered phase is at approximately z = 0.9^ 
indicated by a horizontal line in Fig.0| For V = 1.5, close 
to the charge order phase transition, this value of z = 0.9 
is reached near C ~ 0.33. The distortion close to z = 0.9 
is realized also for other interactions, e.g. deep in the 
ordered phase at y = 3.0, C = 0.24. Since, however, 
NaV205 is probably close to a quantum critical point 
of charge ordering,'^" and because of the value for C ob- 
tained from the band structure calculations, we conclude 
that the system can best be described in the ordered 
phase using V « 1.3 and C « 0.35. 

When calculating the ground-state energy of the sys- 
tem, the question arises how the terms in Eq. ^ con- 
tribute to the total energy, i.e., whether some sort of 
virial theorem holds. In Fig. the behavior of the two 
contributions Ij2b|l and is shown. One can easily see 
that the crossing points of the curves are at the same 
value of z at which the ground-state energy reaches its 
minimum. Fig. |21 Therefore we find that a virial theo- 
rem such as that for the one-electron polaronic stateS)S^ 
is fulfilled in the form 

{Hi)^-\{H,^i), (6) 

with a relative numerical accuracy of better than 10~^. 
For large Coulomb interactions well above the phase 
transition this high accuracy is likely achieved since the 
ion charges depend very weakly on z. Therefore, com- 
pared to He-i and Hi, i?EHM is almost independent of z, 
namely [cLH-e^yim/ dz{z-cai-a) ~ 0.005], and the dependence 
of He-i{z) in Fig. |31is close to a straight line. Then the 
virial theorem follows from the functional form of Hi{z) 
and He-i{z). For smaller values of the interaction, e.g., 

V = 1.5, the dependence of the ion charges and i?EHM 
is considerably larger [dH-EmA/ dz{zm\n) ~ 0.03]. Yet also 
in this case the virial relation is satisfied. The contribu- 
tion of the sum of the lattice terms {Hi + H^-i) to the 
total energy varies between 19% aX V — 1.5 and 30% at 

V = 3.0. 



B. Kink excitations 

So far we have considered the perfect zig-zag charge 
order pattern described by Eq. (Q. This ordering can 
be destroyed by local in-rung excitations where an elec- 
tron hops from the site with minimal energy to that with 
maximal energy, as marked by black and white circles in 
Fig. n The excitation energy of this process is 2V , pro- 
vided the system is totally ordered, that is, mco = l-O- 
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(b) 



FIG. 6: (a) A schematic plot of a non-local kink- like excitation 
in an ordered ladder. The darkness of the circles corresponds 
to the charges on the sites of the ladder, (b) Local kink ex- 
citation with a sharp change of the order parameter. The 
electron states for one of the electrons moving between sites 
u and d are degenerate in this case. 




FIG. 7: Twisted ("Mobius") boundary conditions which pro- 
duce a kink excitation. The numbers label the sites in the 
8x2 cluster from 1 to 16. 




FIG. 8: Ground-state energy per site as a function of the kink 
length L, Eq. 0. 
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FIG. 9: Kink excitation energy as a function of the interaction 
V . It is zero for V < 1.0. The lines are guides to the eye. 



Another type of excitation is the formation of a local 
pair of doubly occupied and empty rungs, which has the 
same energy 2V at mco = l-O- In addition, there are 
nonlocal kinklike excitations, where the order parameter 
smoothly changes along the ladder between two degener- 
ate patterns as shown in Fig. O The nonlocal character 
of the excitation leads to a decrease of the excitation en- 
ergy. Since the lattice is coupled to ion charges by the 
Holstein interaction, the kinks couple to the lattice, too. 

In order to investigate kink excitations we used the 
largest system size available, which is a single ladder 
consisting of eight rungs, and imposed twisted "Mobius" 
boundary conditions'^ as shown in Fig. [7\ The zig-zag 
distortions of Eq. ((SJ are modified to a kink distortion 



ze'^-'^^ tanh 



(Ri — Ro)ea 



L 



(7) 



with the center of the kink Rq located in the middle 
between the rungs, L being its length in units of the 
lattice spacing, and the unit vector in ladder direction. 
We kept z at its previous optimal value Zmin and varied 
L looking for the minimum in the ground-state energy as 
shown in Fig. |S1 



For interactions up to y = 1.0 the ground-state energy 
strongly decreases with increasing kink length L, without 
a minimum, implying that at this weak coupling we have 
no kink excitations in the system. For larger interactions 
we found a clear minimum in the ground-state energy. 
The kink length at V = 1.5 is L w 2.15, and it shrinks 
to L « 1.33 at = 3.0. 

The kink excitation energy is defined as the difference 
between the ground-state energy with twisted boundary 
conditions at the optimal value of L and the ground-state 
energy with periodic boundary conditions and static dis- 
tortions z = Zinin- It is shown in Fig. O In the fully 
ordered state and in the atomic limit where V ^ ta, the 
kink potential energy AE is close to V (see Fig. ^ . The 
actual total energy is considerably smaller since the kinks 
are extended and since at the kink boundary the electron 
energy on the upper (m) and lower (d) legs become de- 
generate as shown in Fig. El and therefore the intrarung 
hopping becomes more likely. The hopping kinetic en- 
ergy of the order of —ta decreases the total energy of the 
system, leading to the result shown in Fig. |51 

For weak interactions, where the kinks are well ex- 
tended and <C 1, our results can be compared with 



6 



a model calculation in a classical 
ladders--- which gives 

1 



model for infinite 
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(8b) 



where Vc is the critical value of the Coulomb interaction 
for the phase transition. To make a connection to our 
results, we estimate Vc from Eqs. H8a|l and (|8bp for the 
distorted lattice at V — 1.5 and C = 0.35 independently 
and compare them. From L and AE we obtain Vc = 1.28 
and Vc = 1.05, respectively. These values are in reason- 
able agreement with each other and consistent with the 
behavior of the order parameter, see Fig. 

At all values of V the kink excitation energy with lat- 
tice coupling is larger than the excitation energy without 
coupling shown as dashed line in Fig.|3 This can be un- 
derstood since at a fixed value oiV > Vc the charge order 
parameter is larger for the distorted lattice. Since the 
ordering is more complete, the kink lengths are smaller 
which increases the excitation energy. Note that without 
lattice coupling we have no parameter L as in Eq. Q for 
the determination of the kink length. 



IV. DYNAMIC PROPERTIES 

In Sec, mil we only considered static properties. How- 
ever, enlightening insight into the physics of the system 
can be extracted from dynamical correlation functions 
showing the spectra of charge and spin excitations. The 
corresponding susceptibilities are given by 



Xc(q,c^) - ke'"*(nq(t)n_q) 



(9a) 
(9b) 



the the Fourier 



where nq(t),n_q and S^{t),S^ 
transforms of the charge and z component of spin densi- 
ties, respectively. 

We calculated the charge susceptibility Eq. on a 
ladder consisting of eight rungs with periodic boundary 
conditions along the ladder. The results are shown in 
Fig. El We define the charge gap Ac as the energy 
at which the lowest lying excitation of the charge sus- 
ceptibility occurs. The corresponding momentum is al- 
ways q = Q. In the disordered phase we have no gap- 
less charge excitation. When increasing the Coulomb 
interaction V, the gap at q = Q decreases as shown 
in Fig. 111! and all other charge excitations become in- 
significant. The charge gap does not vanish exactly for 
C = 0, i.e. without coupling to the lattice, but appears 
to go to zero as N, 



rungs 



oo. When the coupling to 
the lattice is switched on, the gap is exactly zero in the 
ordered phase where the symmetry is broken explicitly. 
The charge gap behaves in a similar way as the order 
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FIG. 10: Charge susceptibility calculated on an 8 x 2 ladder. 
Upper panel; without lattice coupling. Lower panel: With 
lattice coupling {C = 0.35). The wave vector scan consists of 
the two ranges {qa,qb)/T' = (0,0) — > (0, 1) and (1,0) — » (1, 1), 
separated by a tick mark on the horizontal axis. An additional 
broadening of width = 0.1 was used. 



parameter (Fig. |3J), namely, without electron-lattice cou- 
pling it changes smoothly across the phase transition, 
whereas the changes for finite coupling are significantly 
sharper. 

In Fig.^lat V ^ 1.0, a gapless excitation at q Q can 
be seen, which occurs due to a small but finite distortion. 
As already discussed in Sec. IIIII this distortion is finite 
only due to the finite-size effects and should be zero in 
the thermodynamic limit. 

We note that the gap Ac in the charge spectrum is 
different from the one commonly used in DMRG calcula- 
tions, A [Eo{N + 2)-Eo{N)]/2, where Eo{N + 2) and 
Eo{N) arc the ground-state energies for systems consist- 
ing of A^ -f 2 and A^ particles, respectively. Indeed, as a 
function of V, A shows a behavior opposite to Ac, with 
A = in the unbroken phase and A > at large Vi^ 

The spin susceptibility [Eq. (|9bp ] calculated on the 
same system is shown in Fig. 1121 The momentum scan 
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FIG. 11: Charge gap Ac at q = Q of the charge susceptibihty 
as a function of V . Sohd hnes: 8x2 ladder. Dashed lines: 
4x2 ladder. Open symbols: Without lattice coupling. Full 
symbols: With lattice coupling (see text for case V = 1.0). 

consists of two ranges, from {qa,qb) = (0,0) to (0,7r), 
and from (tt, 0) to (7r,7r). In the first range {qa — 0) 
one can clearly see ttie dispersion of an effective one- 
dimensional Heisenberg model, as predicted by pertur- 
bation theory.^^ The main change of the spin suscepti- 
bility as a function of V in this range is a decrease of 
the effective magnetic exchange interaction Jeff with in- 
creasing charge order. It can be extracted from the spin 
dispersion using Jeff = 2aj(0, 7r/2)/7r (Ref. l25|) . and is 
shown in Fig. 1131 According to magnetic susceptibihty 
measurements Jeff in the low-temperature phase is ap- 
proximately 0.8 of the exchange in the disordered phase. 
Our results are in a good qualitative agreement with 
these data in the sense that an increase in charge order 
goes together with a decrease in Jeff-— However, a quan- 
titative comparison cannot be made since in the experi- 
ment the ordering and, correspondingly, Jeff, are traced 
as a function of temperature for given other system pa- 
rameters while we investigate the ordering at T = as 
a function of the extended Hubbard repulsion V. Our 
calculation also agrees well with the analytical results of 
Refs. JJi and where it was shown that the exchange 
rapidly decreases with increasing V. Quantitatively, at 

V = 3,C = our results give Jeff « 0.06 while per- 
turbation theorjiiS*2S predicts Jefi « 0.04. Moreover, for 

V ~ 1.3, C — 0.35, which give a lattice distortion close 
to the value observed experimentally (see Fig. 4), the 
exchange parameter in Fig. 13 is about 67meV which is 
very close to the experimental 60meV observed in the 
inelastic neutron scattering measurements of Ref. 

The second range in the spin spectrum (Fig. I12|l , with 
qa = TT, shows only high energy excitations at small V, 
but again an effective Heisenberg dispersion at large V. 
For small V the gap in the spin spectrum is very close to 
the charge gap, indicating that it is due to charge excita- 
tions. To verify this conjecture, we calculated charge and 
spin susceptibilities in the noninteracting limit V = 0, 
with tb = (isolated rungs). In this case charge and spin 
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FIG. 12: Spin susceptibility calculated on an 8 x 2 ladder. 
Upper panel: without lattice coupling. Lower panel: With 
lattice coupling (C = 0.35). Presentation as in Fig. 1101 

susceptibilities are equal for qa = and the gap is exactly 
the difference between the bonding and the antibonding 
state given by 2ta. Secondly, we analyzed the depen- 
dence of the spin susceptibility on the hopping tf, along 
the ladder in the disordered phase at = 0.5 fFig. I14|l. 
Whereas the dispersion for qa = scales as i^, which 
is clear evidence of the magnetic origin of these excita- 
tions, the difference between the maximal and minimal 
excitation energy for ga = tt scales as ti,. These obser- 
vations show a direct interplay between the spin and the 
dipole-active charge excitations, which is similar to the 
"charged" magnons introduced in Ref. "3^ for interpreta- 
tion of the infrared absorption spectra of NaV205. 

It is interesting to note that the spin spectra in the 
ordered phase appear to possess a mirror symmetry with 
respect to the central tick mark in Fig. 1121 To quantify 
this observation, the dispersions of the low-energy exci- 
tations atV = 3.0 have been depicted in Fig.^]The dis- 
persions for qa — and qa = are indeed very similar. 
Without lattice coupling the dispersion with = tt is 
shifted upwards compared to qa — Q because of the small 
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FIG. 13: Effective magnetic exchange interaction Jeff in lad- 
der direction in units of Jes{V = 0.5) as a function of V, 
extracted from the spin susceptibility Eq. IQbll . The interac- 
tion is shown with (solid line) and without lattice coupling 
(dashed line). 
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FIG. 14: Spin susceptibility in the disordered phase at V = 
0.5 for hopping along the ladder ti, = 0.5 (left) and tt — 0.3 
(right). Momentum scan as in Fig. 1101 



FIG. 15: Spin dispersion in the ordered phase at V = 3.0 
extracted from Fig. 1121 Bottom: with lattice coupling, top: 
without lattice coupling. Direction of momentum scans: Solid 
lines: (0,0) (0,7r), dashed lines: (7r,7r) (7r,0). 



for V = 2.5, C = 0.0 and V = 1.5, C = 0.35, where 
the values of the order parameter are similar, see Fig. |31 
The spin and charge excitations shown in these plots are 
qualitatively the same and the susceptibilities differ only 
slightly on a few points. From these figures we conclude 
that the dynamical spin and charge susceptibilities of the 
system mainly depend on the order parameter but not on 
the way in which the order has been achieved. 



HUBBARD-HOLSTEIN MODEL 



but finite charge gap a,t V = 3.0 (see Fig. [TT|) . With 
lattice coupling and at interactions where no charge gap 
occurs the agreement is even better. This behavior can 
be understood in the following way. In the disordered 
phase where each electron on a rung occupies a molecu- 
lar orbital consisting of two sites, momenta q = (0, tt) and 
q = {tt, 0) are not equivalent (same spin on the two sites 
of the rung, versus opposite spin on two sites of neigh- 
boring rungs). In this phase pure spin excitations with 
qa = TT are not possible since they require different spins 
on different sites within a rung. This could be achieved 
only by exciting another electronic state within the rung, 
which has the energy 2ta. In the totally zig-zag ordered 
state where the electrons are located on one site of the 
rung, these momenta become equivalent. The same holds 
for momenta q = (0, 0) and q = (tt, tt). 

The overall effect of charge ordering on the dynamical 
susceptibilities can best be seen by comparing the plots 



So far we have only considered static distortions of the 
lattice. Although this is a good approximation if the 
dynamic fiuctuations around these equilibrium positions 
are small, quantum phonon effects can play an important 
role, especially in the critical region. In this section, we 
therefore consider the extended Hubbard-Holstein model 
(EHHM) 



H = Hi 



EHM 



E 



Cz 



(10) 



with i/EHM defined in Eq. H2a|) and M being the mass 
of the local oscillators. The operators Zi and pi are the 
coordinate and momentum of the ion on lattice site i, 
and all other quantities are defined in Eq. When 
expressed in phonon creation and annihilation operators, 
it reads (up to a constant) 

H = HEiiM+uJoY.blb,-gY,ibl + hW (11) 
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FIG. 16: (a) Order parameter iriQQ, calculated on a cluster 
with four rungs. Phonon parameters are C = 0.35, At — 0.125, 
and uJo ~ 125 meV. Parameters are chosen as shown in the 
caption, (b) Amco = nico{Nph = 3) - mco(-^ph = 0) for 
two phonon frequencies wo = 60 and 125 meV. Upper two 
curves: Without static distortion. Lower two curves: with 
static distortion. The dotted line marks Amco = 0. 



Here bl (6^) creates (annihilates) a phonon of frequency 
luq (with ?i = 1) at lattice site i, and the phonons are 
locally coupled to the electron density with coupHng 
strength g = C ^JloqI^k,. 

We would like to point out that the nature of the local 
phonon mode in Eq. (|lll) is not specified and could corre- 
spond to one of several phonon modes in the vanadates. 
Clearly, the use of dispcrsionless Einstein phonons ne- 
glects any coupling between lattice distortions of neigh- 
boring sites. However, a coupling of the Holstein type 
is the strongest phonon mode in ^aNiQ^M> It also rep- 
resents the simplest model for electron-phonon interac- 
tions, and has been successfully used to describe the 
physics of other transition metal oxides such as the 
manganitesi^ 

Compared to exact diagonalization of the model de- 
scribed by the Hamiltonian in Eq. , an additional dif- 
ficulty arises in the case of the EHHM since the number 
of phonons is not conserved. Consequently, even for a 
finite number of lattice sites, the Hilbert space contains 
an infinite number of states, and has to be truncated in 
some way in order to apply the Lanczos method. For this 
reason the size of the systems which can be investigated 
is also considerably reduced. We restricted ourselves to a 
lattice with four rungs and chose a subset of the phonon 
states as^^iM 



N ^ (r) 



-1 \iv^:\ 



(12) 



(r) 

where vl denotes the number of phonons at lattice site 
i and |0)ph is the phonon vacuum state. Alternatively, 



the basis states could also be formulated in momentum 
space as in Ref. Is^l 

Now the truncation of the Hilbert space consists of re- 
stricting the total number of phonons in the |r)ph subset 
as 



N 

E 

1=1 



(13) 



leading to {Nph + N ~ l)\/[NphliN - 1)1] allowed phonon 
configurations for N sites. We would like to point out 
that the set of basis states in Eqs. H12|l and consists 
of all possible phonon states with up to A^ph phonons ex- 
cited. In particular, the Hilbert space includes all linear 
combinations of such states . 

Usually, -/Vph is increased until convergence of an ob- 
servable of interest O is achieved. The latter can be mon- 
itored by calculating the relative error |0(A'ph + 1) — 
0(A^ph)|/|0(^ph)|. 

Due to the complexity of the EHHM and the value of 
parameters, it is not possible to include enough phonon 
states to obtain converged results. Nevertheless, as in 
the case of the pure Holstein model^^ it is still possible 
to deduce the tendency of the results as A^ph is increased 
and thereby obtain information about the exact results 
(corresponding to A'ph = oo). 

To reduce the required number of dynamical phonons, 
it is expedient to introduce static distortions Zi as a co- 
ordinate transformation Zi — Zi + Xi so that quantum 
fiuctuations Xi take place around the position Zi. Apply- 
ing this transformation to Eq. Hl()() yields 



H = Hi 



EHM 



2j CZiTLi 



E 



1 K 



(14) 



which in second quantization results in an expression 
analogous to Eq. Hll|l . Note that the first line in the 
above equation is the same as Eq. For the static 

distortions Zi we again use the zig-zag order pattern Q 
and determine the optimal value of z by minimizing the 
ground-state energy in the presence of phonons yielding a 
static distortion Zstat, which is related to Zmin introduced 

in SccHniby Zmin = ^statlA^ph = 0). 

Note that we perform this coordinate transformation 
only because the number of phonons accessible in our 
calculations is very small, and in this case it is better to 
start from a different equilibrium position z = Xstat and 
not from z = 0. If it were possible to use Aph = oo, 
this coordinate transformation would have no influence 
on the physical results and the actual lattice distortions 
would be produced by the dynamical phonons as a coher- 
ent state of oscillators associated with the ions, indepen- 
dent of any initial coordinate transformation. Of course 
the broken symmetry would only occur in the thermody- 
namic limit, while correlations of the phonon positions 
exist already on finite lattices. 
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The effect of dynamical pfionons on the charge order 
parameter is shown in Fig. ^| We did calculations for 
several values of V at phonon frequencies ujq = 60meV 
and ujQ — 125 meV, the two most relevant modes in 
NaV205.^^ The smaller frequency belongs to a collective 
vibration which includes displacements of the vanadium 
and oxygen ions, whereas the larger one corresponds to 
a vibration of the apical oxygen along the z axis. 

From the upper panel of Fig. E| one can easily see that 
for the nontransformed coordinates (circles) the inclusion 
of dynamical phonons with ujq — 125 me V considerably 
increases the charge ordering. Calculations with A^ph = 1 
and 2 (not shown) revealed that the increase is mono- 
tonic in the number of phonon states and we conclude 
that for convergence many more phonon states would be 
necessary. For the distorted lattice (diamonds), the dy- 
namical phonons actually decrease the charge ordering 
for V > 1.0, and the strongest effect occurs in the vicin- 
ity of the phase transition at V = 1.0 and V = 1.5. 
The reason for this decrease is that Zgtat is shifted down- 
wards with increasing number of dynamical phonons. At 
V — 1.0, where a finite Zgtat at iVph = is a finite size 
effect (Sec. IIII.^ . ^stat is reduced to zero for TVph = 3. 
At y > 1.5 the relative change in Zstat is similar to that 
in m^Q fFig. ll6|l . We want to mention that the two solid 
curves in the upper panel of Fig. 1161 give an upper and a 
lower boundary for the actual value of the order parame- 
ter on the four rung lattice, since for A^ph oo results 
for z = and z — Zgtat become equivalent, as discussed 
above. 

In the lower panel of Fig.^|the difference Am^Q of the 
order parameter in the presence of dynamical phonons 
to A^ph = is shown including data for = 60meV. 
For z = (upper two curves) this deviation is positive 
and the effect is always larger for ujq = 125 meV, which 
corresponds to a larger value of g. For z — Zgtat (lower 
two curves) it is negative for V > 1.0. The crossing 
at F = 1.0 is due to the small finite value of Zstat for 
A^ph — 0, i.e., a finite size effect. 

It is interesting to study the order pattern of the dy- 
namically induced distortions. For this purpose we define 
the correlation function 

= E-"'''''*"''^K(-' - (^^))(% - (%)))' (15) 

which measures the zig-zag ordering of the lattice dis- 
tortions, similar to Eq. (0J for the charge densities. For 
the nontransformed coordinates it is depicted in Fig. El 
for ujQ = 125 meV and different numbers of dynamical 
phonons. From this figure it is clear that dynamical 
phonons induce zig-zag lattice distortions which strongly 
increase around the phase transition point. Note that the 
correlation function Cz is not normalized to the interval 
[0, 1] since lattice distortions are not conserved quanti- 
ties, different from, e.g., charges. 

For the transformed coordinates we can calculate the 
dynamically induced zig-zag distortions z^yn directly 
from the expectation values (z^) since the symmetry of 




V 



FIG. 17: Correlation function Cz as defined in Eq. I|15|l with- 
out static distortions for uuo = 125 meV. The number of 
phonons is given in the caption. 

the system is broken explicitly. Results show that the 
sum ztot = ^stat + ^dyn of the static distortion and the 
dynamically induced distortion is always smaller than the 
value Zniin determined in Sec. lIII K\ and this effect is most 
pronounced near the phase transition. For V — 1.5 and 
Nph = 3 we got Zstat — 0.889 and Zdyn = 0.013 yielding 
a total zig-zag distortion of Ztot = 0.902, which is notice- 
ably smaller than Zj^in = 1.001 for Aph = 0. Well above 
the transition point the dynamically induced distortions 
are very small, for instance for V = 3.0 and A^ph = 3 
calculations gave Zdyn = 0.0005, and Zgtat = 1.290 is only 
slightly smaller than Zmin ~ 1.294. 

From this analysis we conclude that using just static 
distortions gives qualitatively correct results but overesti- 
mates the lattice distortions, in particular, in the vicinity 
of the phase transition. The value of the order parameter 
in the full dynamical model with Aph = oo will there- 
fore be somewhat smaller than in Sec. lIII Al as discussed 
above. 



VI. CONCLUSIONS 

In this paper we investigated the influence of lat- 
tice effects on the charge ordering transition, the spec- 
trum of kink excitations and dynamical susceptibilities 
of quarter-filled ladder systems and considered the a'- 
NaV205 compound as an example. For this purpose 
we modified the Hamiltonian of the extended Hubbard 
model by terms which take into account the coupling 
of electrons to lattice distortions. The lattice rigidity k 
and the Holstein coupling constant C used in our model 
were determined by first-principles band-structure cal- 
culations. The physical properties were calculated by 
the exact diagonalization technique. The results for 
the ground-state energy and the order parameter show 



11 



that by including static distortions, the phase transition 
is shifted significantly downward to lower values of the 
nearest-neighbor Coulomb interaction V. The calculated 
displacements of the vanadium ions due to the electron- 
lattice coupling in the charge ordered phase are in good 
agreement with experimental measurements ii We also 
found a virial theorem to be fulfilled to high precision 
for the terms in the Hamiltonian that couple to the lat- 
tice. 

As low-energy excitations of the distorted ground state 
we considered kink excitations, where the charge order 
pattern changes along the ladder between two degener- 
ate configurations. These kinks are long when itIqq is 
small, and become shorter with increasing order. The 
kink lengths and energies at small m^g ^-re comparable 
with those of a classical 0"* modeljiS 

Moreover, we studied the extended Hubbard-Holstein 
model to investigate the effect of dynamical phonons. Re- 
sults showed that they have a strong influence on the 
charge order parameter in the vicinity of the phase tran- 
sition. An analysis of the correlations of the dynami- 
cally induced distortions revealed that phonons indeed 
favor zig-zag lattice distortions. We showed that using 
just static distortions somewhat overestimates the actual 
value of the lattice distortion, but well away from the 



transition point this dynamical effect is very small and 
a description by static distortions gives already accurate 
results. 

In addition to these static properties we also calculated 
the dynamic charge and spin susceptibilities. We showed 
that the main features of these quantities are determined 
by the value of the order parameter and not by the way 
this value is achieved. From the spin susceptibility we ex- 
tracted the effective magnetic exchange interaction along 
the ladder, which exhibits a pronounced decrease with in- 
creasing charge order. The magnitude of this parameter 
taken at V = 1.3, C = 0.35 is in good agreement with 
the experimental inelastic neutron scattering data^ 
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